NASA COOPERATIVE AGREEMENT NCC1-1 12 


STUDY OF THE TRANSITION FLOW REGIME 
USING MONTE CARLO METHODS 


Final Report 


Prepared by 


H. A. Hassan 
Principal Investigator 

Department of Mechanical and Aerospace Engineering 
Campus Box 7910 
North Carolina State University 
Raleigh, North Carolina 27695-7910 



During this period a presentation was made at the 7th AIAA/ASME Joint Thermophysics 
Conference. Another publication, dealing with turbulence at hypersonic Mach numbers, will be 
presented at the AIAA Aerospace Sciences meeting in January 1999. Copies of these publications 
are made part of this final report. In addition, progress made by Professor Bird is enclosed. 

In addition, a list of published research conducted at N. C. State University and supported 
by the cooperative agreement is enclosed. It provides a list of students supported by the 
Cooperative Agreement. 


Publications supported in part by NCC1-1 12 

1. “Grid Generation and Adaptation for the Direct Simulation Monte Carlo Method”, D. P. 
Olynick, J. N. Moss and H. A. Hassan, AIAA Paper 88-2734, June 1988, also Journal of 
Thermophysics and Heat Transfer, Vol. 3, pp. 368-373, October 1989. 

2. “Influence of Afterbodies on AOTV Flows”, D. P. Olynick, J. N. Moss and H. A. Hassan, 
AIAA Paper 89-0331, January 1989. 

3. “Monte Carlo Simulation of Reentry Plasmas”, A. B. Carlson, H. A. Hassan, and J. N. 
Moss, AIAA Paper 89-0683, January 1989. 

4. “Study of Hypersonic Flow Past Sharp Cones”, J. C. Taylor, J. N. Moss and H. A. Hassan, 
AIAA Paper 89-1713, June 1989. 

5 . “Direct Simulation of Reentry Flows with Ionization”, A. B. Carlson and H. A. Hassan, 
AIAA Paper 90-0146, January 1990, also. Journal of Thermophysics and Heat Transfer, Vol. 
6, pp. 400-404, July 1992. 

6. “Monte Carlo Simulation of Vibrational Relaxation in Nitrogen”, D. P. Olynick, J. N. Moss, 
and H. A. Hassan, AIAA Paper 90-1767, June 1990; Also, Journal of Thermophysics and 
Heat Transfer, Vol. 6, pp. 22-26, January 1992. 

7. “Monte Carlo Simulation of Re-Entry Flows Using a Bimodal Vibration Model”, D. P. 
Olynick, J. N. Moss, and H. A. Hassan, Journal of Thermophysics and Heat Transfer, Vol. 
4, pp. 273-277, July 1990. 

8. “Direct Simulation Monte Carlo with Ionization and Radiation”, A. B. Carlson and H. A. 
Hassan, 17th International Symposium on Rarefied Gas Dynamics, Edited by A. E. Beylich, 
VCH, Weinheim, Germany, pp. 763-769, 1991. 


9. “Monte Carlo Simulation of Nonequilibrium Shock Fronts”, D. P. Olynick, J. N. Moss, and 
H. A. Hassan, AIAA Paper 91-1341, June 1991; Also, Journal of Thermophysics and Heat 
Transfer, Vol. 6, No. 4, pp. 626-630, Oct. -Dec. 1992. 

10. “Radiation Modeling with Direct Monte Carlo”, A. B. Carlson, and H. A. Hassan, AIAA 
Paper 91-1409, June 1991; Also, Journal of Thermophysics and Heat Transfer, Vol. 6, No. 
4, pp. 631-636, Oct.-Dec. 1992. 

1 1 . “Radiation Transport Around Axisymmetric Blunt Body Vehicles Using a Modified 
Differential Approximation”, Lin C. Hartung and H. A. Hassan, AIAA Paper 92-0119, 



January 1992; Also, Journal of Thermophysics and Heat Transfer, Vol. 7, No. 2, pp. 220- 
227, April-June 1993. 

12. “Monte Carlo Simulation of Reentry Flows with Ionization”, Jeff C. Taylor, Ann B. Carlson, 
and H. A. Hassan, AIAA Paper 92-0493, January 1992. 

13. “Monte Carlo Simulation of Entry in the Martian Atmosphere”, David B. Hash and H. A. 
Hassan, AIAA Paper 92-0494, January 1992; Also, Journal of Thermophysics and Heat 
Transfer, Vol. 7, No. 2, pp. 228-232, April-June, 1993. 

14. “Direct Simulation with Vibration-Dissociation Coupling”, D. B. Hash and H. A. Hassan, 
AIAA Paper 92-2875, July 1992; also Journal of Thermophysics and Heat Transfer, Vol. 7, 
No. 4, pp. 680-686, Oct.-Dec. 1993. 

15. “A New Two-Temperature Dissociation Model for Reacting Flows”, D. P. Olynick and H. A. 
Hassan, AIAA Paper 92-2943, July 1992; also Journal of Thermophysics and Heat Transfer, 
Vol. 7, No. 4, pp. 687-696, Oct.-Dec. 1993. 

16. “Direct Simulation of Diatomic Gases Using the Generalized Hard Sphere Model”, D. B. 

Hash and H. A. Hassan, AIAA Paper 93-0730, January, 1993. 

17. “A Generalized Hard-Sphere Model for Monte Carlo Simulation”, H. A. Hassan and D. B. 
Hash, The Physics of Fluids, Vol. 5, No. 3, pp. 738-744, March 1993. 

18. “Monte Carlo Simulation of Radiating Re-Entry Flows”, Jeff C. Taylor, Ann B. Carlson, and 
H. A. Hassan, AIAA Paper 93-2809, July 1993; Also, Journal of Thermophysics and Heat 
Transfer, Vol. 8, No. 3, pp. 478-485, July-September 1994. 

19. “Comparisons Between DSMC and Navier-Stokes Equations for Re-Entry Flows”, D. R. 
Olynick, J. C. Taylor, and H. A. Hassan, AIAA Paper 93-2810, July 1993, Also, Journal of 
Thermophysics and Heat Transfer, Vol. 8, No. 2, pp. 251-258, April-June 1994. 

20. “Monte Carlo Simulation Using Attractive-Repulsive Potential”, D. B. Hash and H. A. 
Hassan, in Rarefied Gas Dynamics: Theory and Simulations, Progress in Astronautics and 
Aeronautics, edited by B. D. Shizzal and D. P. Weaver, Vol. 159, pp. 284-293, 1994. 

21. “Effects of Chemistry and Rarefaction on Blunt Body Wake Structure”, V. K. Dogra, R. G. 
Wilmoth, J. N. Moss, J. C. Taylor, and H. A. Hassan, AIAA Paper 94-0352, January 1994; 
Also, AIAA Journal, Vol. 33, No. 3, pp. 463-469, March 1995. 

22. “Rates of Thermal Relaxation in Direct Simulation Monte Carlo Methods”, B. L. Haas, D. B. 
Hash, G. A. Bird, F. E. Lumpkin, III and H. A. Hassan, Physics of Fluids, Vol. 6, No. 6, 
pp. 2191-2201, June 1994. 

23. “Upwind Methods for Ionized Flows”, J. C. Taylor, D. R. Olynick, and H. A. Hassan, 
AIAA Paper 94-1957, June 1994. 

24. “Rarefaction Effects on Blunt Body Wake Structure for Earth Entry Conditions”, V. K. 

Dogra, J. N. Moss, R. G. Wilmoth, J. C. Taylor, and H. A. Hassan, AIAA Paper 94-2016, 
June 1994; Also, Journal of Thermophysics and Heat Transfer, Vol. 9, No. 3, July 
September 1995, pp. 464-470. 

25. “Parallel DSMC Solution of Three-Dimensional Flow Over a Finite Flat Plate”, R. P. Nance, 
R. G. Wilmoth, B. Moon, H. A. Hassan, and J. Saltz, AIAA Paper 94-2019, June 1994; 
Also, Journal of Thermophysics and Heat Transfer, Vol. 9, No. 3, July-September 1995, pp. 
471-477. 



26. “Direct Simulation of Diatomic Gases Using the Generalized Hard Sphere Model”, D. B. 
Hash, J. N. Moss, and H. A. Hassan, Journal of Thermophysics and Heat Transfer, Vol. 8, 
No. 4, Oct.-Dec. 1994, pp. 758-764. 

27. “A DSMC Navier-Stokes Hybrid Solver”, David B. Hash and H. A. Hassan, AIAA Paper 
95-0410, January 1995. 

28. ‘The GHS Interaction Model for Strong Attractive Potentials”, J. A. Kune, D. B. Hash, and 
H. A. Hassan, The Physics of Fluids, Vol. 7, No. 5, pp. 1173-1175, May 1995. 

29. “A Radiation Absorption Model for the DSMC Method”, J. C. Taylor, H. A. Hassan, and L. 
Hartung Chambers, in Rarefied Gas Dynamics 19, J. K. Harvey and R. G. Lord, Editors, 
Vol. 2, 1995, pp. 829-835. 

30. “Rarefaction Effects on Blunt-Body Wake Structure”, V. K. Dogra, J. N. Moss, J. C. 

Taylor, D. B. Hash, and H. A. Hassan, in Rarefied Gas Dynamics 19, J. K. Harvey and R. 

G. Lord, Editors, Vol. 2, 1995, pp. 1154-1160. 

31. “Direct Simulation of Shock Front Radiation in Air”, A. K. Berghausen, J. C. Taylor, and H. 
A. Hassan, AIAA Paper 95-2051, June 1995; Also, Journal of Thermophysics and Heat 
Transfer, Vol. 10, No. 3, July-September, 1996, pp. 413-418. 

32. “A Decoupled DSMC/Navier-Stokes Analysis of a Transitional Flow Experiment”, D. B . 
Hash, and H. A. Hassan, AIAA Paper 96-0353, January 1996. 

33. “A Comparison of Grid-Definition Schemes for DSMC”, R. P. Nance, R. G. Wilmoth, and 

H. A. Hassan, AIAA Paper 96-0604, January 1996; also. Journal of Thermophysics and 
Heat Transfer, Vol. 11, No. 2, April-June, 1997, pp. 296-303. 

34. “Assessment of Schemes for Coupling Monte Carlo and Navier-Stokes Solution Methods,” 
D. B. Hash and H. A. Hassan, Journal of Thermophysics and Heat Transfer, Vol. 10, No. 2, 
April-June 1996, pp. 242-249. 

35. “A Novel Approach for Calculating Equilibrium Radiating Flows,” F. Brauns and H. A. 
Hassan, AIAA Paper 96-1887, June 1996; also. Journal of Thermophysics and Heat 
Transfer, Vol. 11, No. 1, January-March, 1997, pp. 52-58. 

36. “Role of Boundary Conditions in Monte Carlo Simulation of MEMS Devices,” R. P. Nance, 
D. B. Hash and H. A. Hassan, AIAA Paper 97-0375, January 1997; Also, Journal of 
Thermophysics and Heat Transfer, Vol. 12, No. 3, July-September 1998, pp. 447-449. 

37. ‘Two-Dimensional Coupling Issues of Hybrid DSMC/Navier-Stokes Solvers,” D. B. Hash 
and H. A. Hassan, AIAA Paper 97-2507, June 1997. 

38. ‘Transition and Turbulence Modeling for Blunt-Body Wake Flows,” R. P. Nance, T. J. 
Horvath, and H. A. Hassan, AIAA Paper 97-2570, June 1997. 

39. “Solution of Transitional Wake Flows at Mach 10,” R. P. Nance, B. R. Hollis, T. J. 
Horvath, and H. A. Hassan, AIAA Paper 98-2939, June 1998. 

40. ‘Turbulence Modeling of Shock-Dominated Flows with a k-C, Formulation,” R. P. Nance and 
H. A. Hassan, AIAA Paper 99-0153, January 1999. 



Subcontract 87-0334-10 to GAB Consulting Pty. Ltd. 

NASA Grant NAG-1-112 to North Carolina State University 

FINAL REPORT: 1 November 1997 to 31 October 1998 

G.A. Bird 

1. New DSMC Procedures 

(a) Transient adaptive sub-cells 

Several recent applications of the DSMC method have shown that unsteady two- 
dimensional flows can now be studied without resorting to ensemble averaging. For 
Knudsen numbers well under 0.01, this requires the flow output to be based on comput- 
ational cells that are large in comparison with the mean free path. However, a basic 
requirement of the DSMC method is that the mean spacing between collision partners 
should be small in comparison with the local mean free path. The two requirements can 
be reconciled by the specification of a large number of sub-cells in each cell. Sub-cells are 
required only bv the collision routine and the latest DSMC programs employ ‘transient’ 
sub-cells in the form of a rectangular grid that is temporarily overlaid on each 
computational cell when collisions are to be calculated. The number of divisions in each 
cell is adapted to the number of simulated molecules in the cell at the time such that 
there is only one or two molecules in each grid location. The molecules are indexed to 
this grid in the usual way and the selection process effectively yields nearest-neighbor 
collisions in a computationally efficient manner. 

These transient adaptive sub-cells have been implemented in Version 3 of the DS2G 
program and have proved to be highly advantageous. Not only is less computation time 
and memory required for a given mean separation distance of the collision partners, but 
the new procedures lead to far smaller separation distances than the smallest possible 
distances with the traditional sub-cells. A further advantage is that the user is not only 
relieved of the responsibility for choosing the sub-cell size, but there is no longer a need 
to set different cell and sub cell sizes over complex flowfields. Instead, the smallest 
possible spacing of the collision partners is automatically achieved in all regions of the 
flowfield. 

In order to check the efficacy of the transient sub-cells in the DS2G program, diagnostic 
output was added to the TECPLOT file for the flowfield properties. This displayed the 
ratio of the local mean separation between collision partners to the local mean free path. 
The ratio of the time step to the local mean collision time was also displayed and this 
highlighted a problem that has been well known, particularly for hypersonic flows with 
large variations in density. If the time step is smaller than the mean collision time in the 
regions of highest density and/or temperature, it is much smaller than the values that 
are needed in the other regions of the flowfield and the calculation is inefficient. This has 
led to a tendency to set an excessively high time step. The general guideline is to set the 
time step such that the number of collisions in the step is about one quarter the number 
of simulated molecules, but the diagnostic output showed that this is much too large for 
problems with large variations in density or temperature. Cell-dependent time steps 
have been successfully employed by some workers, but these lead to problems with 
efficient programs that do not calculate trajectory intersections with cell boundaries. 
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(b) Molecule and cell dependent time steps 


New procedures have been introduced that make the time step vary with every molecule 
as well as with every cell . To this end, separate time variables are associated with every 
molecule and every cell. The procedures keep all these variables concurrent with the 
overall time variable. 

The overall time variable is advanced in very small time steps, but only a small fraction 
of the molecules are moved and collisions are calculated in only a small fraction of the 
cells at any one time step. The mean collision time is kept for each cell and collisions in a 
cell appropriate to a time interval of 0.3xthe local mean collision time are calculated 
whenever the cell time variable falls 0.15xthe local mean collision time behind the overall 
flow time. Similarly, a molecule will generally be moved through 0.3xthe mean collision 
time of its current cell when its time variable falls 0.15xthis local mean collision time 
behind the overall flow time. However, there may be other constraints on the length of 
the molecule time step based. These include the well-established limit on the time step 
near the axis when “two-dimensional logic” is employed in axially symmetric calculations. 
Also, the logic can often be simplified if there is a molecule velocity dependent restriction 
that limits the distance that a molecule can move in a single step. The step through 
which the overall time is advanced is kept under 0.15xthe minimum mean collision time 
in any cell. The values of these multiples will be optimized as experience is gained with 
the new procedures. 

Because of the very small overall time step and the small number of cells for which 
collisions are calculated at each step, the implementation of the traditional indexing 
routine at each step would involve a prohibitive computational effort. This problem has 
been solved through the introduction of “bi-directional” indexing. The traditional index 
array is comprised of the molecule numbers in order of the cells. If the address of the 
corresponding entry in this array is stored for each molecule, the indexing may be kept 
up to date continuously as molecules move from cell to cell. The index array must also be 
extended with sufficient free space between cells to cope with the fluctuations in the 
number of molecules in individual cells. Some general rearrangement takes place from 
time to time to allow for particularly large fluctuations or systematic changes in the 
distribution of molecules. An additional but minor change is that the entering molecule 
flux must be associated with the cells and the entry of new molecules must occur during 
the collision routine rather than during the move routine. In many other procedures, the 
ability to assign different times to molecules leads to a simplification of the program. 

These new procedures have been implemented in a “test” version of the DS2G program. 
The computational speed of the first demonstration case increased by a factor of more 
than three. This was for a flow with a stream Mach number of two and a realistic 
hypersonic test case showed a speed gain factor between four and five. Other test cases, 
such as uniform flows, were calculated and there was no sign of undesirable side effects. 

The new procedures are not only more efficient but, because of the automatic setting of 
the time step, it is no longer possible for this to be incorrectly specified. The transient- 
adaptive sub-cells effectively prevent the incorrect setting of cell size and a poor 
calculation can be made only if the total number of simulated molecule is too small. This 
could cause the mean separation of collision partners to exceed the local free path but, 
since this is continuously checked, the program can be made to stop if a bad calculation is 
attempted. These new procedures enable DSMC programs to be made foolproof. 
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2. The DS3W and DS2A Programs 


The introduction of the transient-adaptive sub-cells permits the use of the simplest 
possible cell structure - a uniform rectangular grid. However, there are difficulties when 
a thin body causes sections of a single rectangular element to be on the opposite sides of 
the body. In addition, there are demonstrable advantages at low Knudsen numbers in 
having surface cells that conform to the geometry. The solution is to wrap cells layers of 
uniform thickness around all surfaces, with the total depth of this surface cell layer set to 
be larger than the basic rectangular cell size. This scheme requires the calculation of the 
distance of a molecule from the surface when it is within the surface cell layer. Fast 
scaling can then be used to determine the layer within a cell group. Also, for steady 
flows, it is possible to implement the “DTM-FNUM” scaling that has proved to be very 
advantageous for entry type calculations. 

For three-dimensional flows over surfaces that are defined by triangular elements, 
surface cell layer groups are defined by their distances above the triangle, above edges 
that are sufficiently convex to the flow, and above vertices that are sufficiently convex to 
the flow. Depending on the degree of convexity, the edge-based layer may have to be 
further divided on the basis of the elevation angle and vertex-based layers may have to 
be divided in both elevation and azimuth. The solid geometry is not trivial, but 
expressions have been derived for the definition of all these types of surface defined cells. 
The junctions between cell groups can be define analytically, but it is computationally 
more efficient to allow some degree of ambiguity and, when more than one cell is valid, 
the closest one is chosen. 

This geometry scheme is being implemented in a program DS3W for the flow past three- 
dimensional bodies. The program is in structured FORTAN 90 and is being written such 
that it can be progressively debugged, subroutine by subroutine, while it is being written. 
While the geometry and cell definition routines have been completed, the introduction of 
the variable time step procedures has created an urgent need for a new general program 
for two-dimensional and axially symmetric flows. Work on DS3W has therefore been 
suspended until program DS2A has been completed. 

Program DS2A employs the two-dimensional version of the geometry in program DS3W, 
but is more general in that there can be multiple surfaces. It can therefore be applied to 
internal as well as external flows. The user will specify the geometry, the flow conditions 
and the number of megabytes of memory that are available. The computational 
parameters will be set automatically and, if the available memory does not allow a 
sufficient number of simulated molecules for a good calculation, the program will stop 
with an appropriate warning. To date, the geometry and cell definition subroutines have 
been completed for this program also. 

3. Papers Presented 

The following conference papers were presented during this period: 

Knudsen and Mach Number Effects on the Development of Wake Instabilities, ALAA 36 th 
Aerospace Sciences Meeting, Reno, January 1998. 

Physics by Default - the Art of Phenomenological Modeling, 5 th International Symposium 
on the Mathematical Aspects of Fluid and Plasma Dynamics, Maui, June 1998. 

Shock Wave Structure in a Gas with Transverse Gradients, 21 st International Symposium 
on Rarefied Gas Dynamics, Marseille, July 1998. 
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4. Miscellaneous Applications and Program Modifications 

Check calculations for Mars entry flows indicated that the procedures in the DS2G 
program for vibrational relaxation rates were less than satisfactory for carbon dioxide. 
The variable relaxation rates had been linked to the macroscopic temperatures in the 
cells. This was satisfactory with gases such as nitrogen and oxygen, but in the case of 
carbon dioxide the large number of internal degrees of freedom produced much larger 
differences between the translational and overall temperatures. This led to anomalous 
relaxation rates and the procedures were modified to make all relaxation processes 
dependent on the energy in the individual collisions. This was a reversion to earlier 
procedures in which dependence on macroscopic temperatures had been avoided on the 
grounds that it was physically unrealistic. The vibrational relaxation behavior of carbon 
dioxide was then consistent with expectations. 

The transient sub-cells in the general program DS2G made it possible for this program to 
deal with unsteady flows without resort to ensemble averaging. The necessary changes 
to the sampling procedures were incorporated into version 3.1. A number of CFD and 
DSMC studies have recently been made of hypersonic biconic flows that include an 
extensive region of separated flow. While these studies have assumed steady flow, the 
flow conditions bordered on those at which unsteady flow would be expected. The 
unsteady DS2G calculations have shown that the unsteady flow development phase is 
more protracted than had been thought and may involve secondary vortices and some 
oscillation of the shock location. The 25°-65° biconic, M= 10 test case is currently being 
calculated by Version 3.1 of DS2G at a Reynolds number approaching 50,000. The single- 
vortex separated region that was characteristic of the lower Reynolds number cases has 
been replaced by a multi-vortex system with the first secondary vortex appearing on the 
25° surface well upstream of the corner. This calculation on a 450 MHz Pentium II 
system with 512 Mb memory will last 10 weeks. The introduction of the new DS2A 
program is expected to cut the execution time for this type of flow by a factor of ten. It 
should then be possible to extend these calculations into the regime where the separated 
region becomes permanently chaotic, although the assumption of axially symmetric flow 
would then be dubious. 

The demonstration of a transition to a permanently chaotic state in axially symmetric 
Taylor-Couette flow has been the most important result from the flow studies under this 
series of subcontracts. Stefanov, Roussinov, Cercignani, Giurin and Struckmeier have 
recently shown (21“ RGD Symposium, Marseille, 1998) that a similar transition occurs in 
two-dimensional Rayleigh-Benard flows. Both DSMC and Navier-Stokes solvers were 
employed in this study and a similar transition to a permanently chaotic state occurred 
for both the particle and continuum approaches. They characterized this behavior as a 
“strange attractor” and state unequivocally that it represents the first stage of the 
transition to turbulence in the flow. To be physically meaningful, transition calculations 
must be three-dimensional. The N-S and DSMC computations were of similar magnitude 
in this case but, since N-S calculations do not scale linearly with the number of grid 
points, DSMC could well be faster when three-dimensional calculations are made. Also, 
because the ratio of the minimum vortex size to the mean free path is expected to be 
inversely proportional to the Mach number, the Navier-Stokes equations may not be 
valid for transition in the case of hypersonic flow. It will be ten to fifteen years before 
three-dimensional DSMC studies of the initiation of turbulence become possible on 
personal computers. These calculations are possible on the largest of the contemporary 
parallel machines, but the very long runs that would be required would probably not be 
cost effective. The calculations will become feasible within the next few years. 
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Solution of Transitional Wake Flows at Mach 10 
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Modeling enhancements are investigated for transition and turbulence in hypersonic 
blunt-body flows. The current approach combines the (Enstrophy) compressible tur- 
bulence model with a new transition mechanism appropriate for free shear layers, such 
as those encountered in blunt-body wakes. Prediction of transition onset is provided 
by a minimum shear-stress criterion applied along the sting. This method is applied to 
several perfect-gas Mach 10 flows over 70-deg blunted cones. Improved prediction of the 
peak sting heating rates is indicated when compared to results obtained using streamwise 
instability modes. Grid refinement is also investigated, and is found to provide addi- 
tional improvement in the agreement between the present approach and experimental 
data. The favorable comparisons are a further indication that transition to turbulence is 
indeed present in the blunt-body wake flows examined. 


Nomenclature 

k Turbulent kinetic energy 

p Pressure 

Q Reynolds- or time-averaged value of Q 

Q Favre-averaged value of Q 

q Heat-transfer rate 

Nose radius 

Re Reynolds number 

s Linear surface distance 

T Temperature 

U Velocity magnitude 

Ui Velocity vector 

x Distance measured along sting support 

Boundary- or shear-layer thickness: 
k - £ model constant 

Boundary-layer displacement thickness 
Sij Kronecker delta 

r Intermittency 

A Transition-extent parameter 

P Dynamic viscosity 

^ Kinematic viscosity 

p Density 

t Characteristic time scale 
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Tij Reynolds stress tensor 

u Transitional frequency 

uji Vorticity vector 

£ Enstrophy 

Subscripts 

1 Value at edge of shear layer 

D Value based on cone diameter 

e Value at edge of boundary layer 

T Total value 

t Turbulent value 

tr Transitional value 

oc Freestream 

Superscripts 

/ Laminar (non-turbulent) value 

t Turbulent value 

♦ Dimensional transitional quantity 

Introduction 

Current and proposed planetary exploration mis- 
sions. such as the recent Mars Pathfinder project, have 
spurred renewed interest in the physics of blunt-body 
wake flowfields. Accurate characterization of the near- 
wake environment is important for the design of entry 
configurations, since the nature of the wake closure 
typically places constraints on payload size and loca- 
tion. Recent activity in this discipline has included a 
number of experiments on blunt-cone models, such as 
the work carried out by several teams of researchers 1 * ^ 4 
in support of AGARD Working Group 18. These WG 
18 tests have been conducted in severed hyperveloc- 
ity facilities with the objective of characterizing the 
fluid dynamic phenomena present in the wake region 
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Fig. 1 Schematic of hypersonic blunt-body wake flow 


for flows exhibiting real-gas behavior (Fig. 1). The 
tests were conducted at nominally identical test con- 
ditions to help assess and quantify facility-to-facility 
performance. Initial comparisons of measurement with 
prediction indicated that the measured peak heat- 
ing rates along the model sting support were two to 
three times greater than those predicted by laminar 
Navier-Stokes computations. A— variety of explana- 
tions were proposed (such as flow establishment, real- 
gas phenomena, and wake rarefaction) to explain this 
disparity. Subsequent perfect-gas tests performed by 
Horvath and coworkers 5 ’ 8 on the AGARD WG 18 ge- 
ometry, as well as studies by Hollis and Perkins 7 using 
a similar configuration, suggested instead that wake 
shear layer transition was responsible for the higher- 
than-anticipated heating levels at reattachment. 

Accurate predictions of such high-speed transitional 
flows are important for the design of hypersonic ve- 
hicles. However, there has been little effort to date 
towards development of computational tools capable 
of predicting transitional and turbulent flow in blunt- 
body wakes. 

Previous efforts by Nance, Horvath and Hassan 3 to 
predict transitional blunt-body flow at a freestream 
Mach number of 6 yielded improved agreement with 
experimental heat-flux data obtained by Horvath. 
McGinlev and Hannemann. 5 In the present work, we 


use the same general approach as in Ref. 8, but the 
mechanism responsible for transition is different. The 
earlier work utilized streamwise instabilities, namely 
Tollmien-Schlichting and Mack modes. Currently, 
however, we consider transition as a result of free- 
shear-layer instability. Two configurations shall be 
considered here, both of which were tested in the 31” 
Mach 10 Air Tunnel at NASA Langley Research Cen- 
ter. These geometries, shown in Figures 2 and 3, 
correspond to those tested by Horvath and Hanne- 
mann 6 and Hollis and Perkins, 7 respectively. 

Modeling Approach 

This work is based on solution of the Favre- averaged 
Navier-Stokes equations in conjunction with a two- 
equation model modified for use in transitional flow- 
fields. In this context, we seek to determine the in- 
fluence of the transitional or turbulent fluctuations 
on the mean flow variables, instead of predicting the 
fluctuations themselves. Closure is obtained by speci- 
fying equations for the modeled variables, as well as 
expressions for disturbance-related quantities in the 
non-turbulent portions of the flow. 

Transport Equations for Modeled Quantities 

Fluctuation velocity and time scales are provided in 
this work by the compressible k-C, model first proposed 
by Robinson, Harris and Hassan. 9 This approach is 
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R b =7.62 cm: RyR b =0.5: R7R b =0.05 
R/R b =0.083: L/R b =6: R JR b =0.25 



Fig. 2 70-deg blunt cone and sting geometry employed by Horvath and Hannemann 


very similar as that used in Ref. 8, with minor modifi- 
cations to some of the terms and constants appearing 
in the model. The k-C equation set used in this work is 
listed below, and model constants are listed in Table 1. 
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Specific alterations to the k-Q equation set include 
reformulation of the density-gradient timescale and re- 
duction of the pressure- work constant C k . These mod- 
ifications were made to preserve Galilean invariance 
while maintaining good agreement with experimental 
data when applied to transonic airfoils. Note that the 
Reynolds stress tensor is modeled using the Eoussinesq 
approximation. Preliminary investigations of these 
flowfields employing a partial differential stress model, 
wherein the Reynolds normal stresses are predicted 
using model transport equations, yielded no improve- 
ment in predictions when compared to the Boussinesq 
approach. 
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R b =2.54 cm: R n /R b =0.5; R./R b =0.05; R/R b =0.60 
R/R^O.406: L s /R b =5.25 



Fig. 3 70-deg blunt cone and sting geometry employed by Hollis and Perkins 


Table 1 k-£ model constants 


Constant 

Value 

Q 

0.09 

a 3 

0.35 

04 

0.42 

05 

2.37 

06 

0.10 

07 

1.4Q-- 


0.065 

<?r 

0.07 


65.0 

l/CTk 

1.80 

l/*< 

1.46 

5 

0.10 

Ci 

0.60 

c k 

1.00 

. 

2.10 


Instability Mechanism 

An advantage of the k - £ formulation is that the 
model terms are derived without explicitly declaring 
the nature of the fluctuation. Hence, the model may 
be used to predict the entire flowfieid, including transi- 
tional and fully turbulent regions, provided an appro- 
priate stress-strain relation is employed. Such a rela- 
tionship is obtained in part by defining the total eddy 


viscosity in terms of a fully turbulent contribution and 
a contribution due to non-turbulent fluctuations: 

Mr = (1 — T) fi t r 4- Ffi t (3) 

Eq. 3 can be rewritten as a combination of terms 
of non-turbulent and fully turbulent timescales as fol- 
lows: 

Mr = C^pkr^ (4) 

r u = (l-r)r' +IY‘ (5) 

where the fully turbulent timescale is given as 


and the non-turbulent timescale is dictated by the 
mechanism responsible for transition. In Ref. 8, 
non-turbulent disturbances were considered through a 
model accounting for Mack’s 10 first and second modes. 
Hence, the non-turbulent characteristic timescale was 
originally defined as 

rf, = n + r 2 (7) 

where the two timescales are defined in terms of first- 
and second-mode characteristic frequencies, and the 
frequencies are correlated as 
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Fig. 4 Base plane and sting heating results ob- 
tained using earlier model-Mach 6, Reo = 2 x 10 s 
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where U v , the phase velocity, is taken to be 0.94 times 
the edge velocity U t . As illustrated in Ref. 8, this 
mechanism yielded good agreement between compu- 
tation and the experimental heating data 5 for a repre- 
sentative Mach 6 case, also using the AGARD WG 18 
geometry. JYeestream conditions for the Mach 6 cases 
discussed here are listed in Table 2. The Mach 6 result 

Table 2 Freest ream conditions for representative 
Mach 6 cases 


Reo 

Poo,kg/m J 

Uoo, m/s 

Too, K 

2 x 10 ti 

6.45 x 10" 2 

943.8 

62.7 

0.5 x 10 a 

1.68 x lO" 2 

931 

63.3 


shown in Fig. 4, which was obtained using the original 
instability model, demonstrates that this approach is 
capable of providing reasonable predictions of the lo- 
cation and magnitude of the peak heat-transfer rate 
adong the sting for this case. The vertical arrow in 
the figure denotes the location of transition onset pre- 
dicted by the computation. When this model is ap- 
plied to representative Mach 10 flowfields, however, 
the agreement is considerably poorer, as illustrated in 
Figures 5 and 6. The degradation of agreement with 
increasing Mach and decreasing Reynolds number is a 
sign of an improperly modeled transition process. 

A correct transition mechanism should yield good 
results for both Mach 6 and Mach 10 experiments over 
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Fig. 5 Base plane and sting heating results ob- 
tained using earlier model-Mach 10, Reo = 1 x 10 6 



Fig. 6 Base plane and sting heating re- 
sults obtained using earlier model-Mach 10, 
Rep = 2.8 x 10 a 

a range of Reynolds numbers. Because the mechanism 
used in Ref. 8 was inadequate for the Mach 10 flow- 
fields, the new model focuses instead on transition in 
the free shear layer emanating from the shoulder of 
the cone, shown in Fig. 1. Such free shear layers are 
always unstable. The most amplified frequency, which 
is required for the present formulation, is obtained 
from linear stability theory predictions for compress- 
ible mixing layers. 11 For all the Reynolds numbers, 
Mach numbers, and and velocity ratios considered in 
Ref. 11, the nondimensionad value of the most am- 
plified frequency was found to be relatively invariant, 
ailways falling in the range u = 0.05 to lj = 0.15. This 
nondimensionad frequency is related to the dimensional 
value by 

n 
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Table 3 Freestream conditions and grid dimensions for Mach 10 cases 


Case 

Geometry 

Hen 


LU Q , m/s 

K 

Grid 1 

Grid 2 

1 

Horvath and Hannemann 

L :< 10° 

1.613 x ur 7 

1426 

31.27 

139 x 75 

377 x 149 

■> 

Horvath and Hannemann 

2.30 x 10* 

4.471 x LIT 3 

1415 

53.15 

189 x 75 

377 x 149 

3 

Hollis and Perkins 

3.15 x 10 6 

1.646 x 1(T 2 

1425 

51.48 

125 x 90 

249 x 179 

4 

Hollis and Perkins 

3.23 x 10' 

4.51 x 10"' 

1416 

53.31 

125 x 90 

249 x 179 


a/ =u>u;i <r 


The intermittency correlation due to Dhawan and 
(11) Narasimha iJ is retained from the previous work: 


where (J\ is the velocity of the higher- velocity stream, 
and 5* is proportional to the thickness of the mixing 
layer: 


<5* = (12) 

In the above, z* represents the dimensional distance 
from the start of the mixing layer, which in this case 
is taken to correspond to distance along the sting sup- 
port. The dimensional frequency can be used to define 
the non-turbulent characteristic timescale: 



(13) 


For the present work, the nondimensional frequency 
appearing in Eq. 11 is taken as 0.10, and the model 
constant d is chosen as 5.0. 

The dissipation term in the k equation is also mod- 
ified to account for transitional effects by computing 
the representative decay time for turbulent kinetic en- 
ergy as 


where 


1 - 1 - r L 

T k ~ Tfc n 



S — S{j Sij 


(14) 


and the model constant C is chosen as 0.01. While 
a more rigorous approach would correlate u with pa- 
rameters such as Mach number and Reynolds num- 
ber, there is insufficient data for performing such a 
correlation at present. The values listed above were 
calibrated through application to the Mach 6, Reo = 
2 x 10 6 blunt-body flow previously considered in Ref. 8. 


Transition Prediction and Intermittency 
Calculation 

As in Ref. 8, we choose the transition onset loca- 
tion as the location of minimum shear stress along 
the sting. The minimum heat-flux criterion proposed 
earlier is not considered here, since it was found to 
yield poorer agreement with available data. A crite- 
rion based on location of flow reattachment along the 
sting was also considered, but was discarded because 
it also yielded degraded comparisons. 


r = 1 — exp (— 0.412£ 2 ) (15) 

However, the intermittency determination in Ref. 8 
utilized a distance parameter £ based on linear surface 
distance starting from the forward stagnation point. 
In keeping with the notion of shear-layer instability, 
the parameter £ is now based on distance along the 
sting. That is, 


, max (x — it,0) 

A (16) 

where 

Re> = 9.0Re| ( 75 (17) 

and the Reynolds numbers appearing above are based 
on conditions at the edge of the layer. In Eq. 16, x t 
is the location of transition onset, and is determined 
as part of the solution to correspond to the location of 
minimum skin friction along the sting. 

Numerical Method 

The modeling proposed here has been incorpo- 
rated into Olynick’s 13 two-dimensional/axisymmetric 
implicit solver for hypersonic flows. This algorithm 
solves the governing equations for 5-species air in ther- 
mochemical nonequilibrium. The solver uses Roe’s 
flux-difference splitting 14 for the inviscid flux, ex- 
tended to higher order using MUSCL variable extrap- 
olation 15 with a minmod slope limiter. Time integra- 
tion is accomplished using the diagonal implicit vari- 
ant of Yoon and Jameson’s Lower-Upper Symmetric 
Gauss-Seidel method, 16 which only requires the inver- 
sion of diagonal matrices. This property is attractive 
for nonequilibrium flows, where a large number of par- 
tial differential equations must be solved. 

Since the current flowfields under investigation are 
perfect-gas, vibrational relaxation and chemical re- 
actions are disabled to obtain the results in this 
study. Additionally, the high-temperature transport- 
property calculations originally used in the code were 
replaced with Sutherland’s law for viscosity and a con- 
stant laminar Prandtl number of 0.72. Closure for the 
Favre-averaged energy equation is accomplished using 
a constant turbulent Prandtl number of 0.89. 

To promote physically realizable solutions, the com- 
puted values of k and Q are required to remain pos- 
itive throughout the computational domain. How- 
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a) Forebody and shoulder 



b) Base plane and sting 

Fig* 7 Influence of grid refinement for case 1— 
Mach 10, Rep = 1 x 10® 

ever, no artificial limitation on turbulent kinetic en- 
ergy production was imposed in the numerical method. 
Transitional solutions are started by running the solver 
in a laminar-flow mode for about 100 iterations. 

Results 

As mentioned earlier, two blunt-cone configurations 
were considered in this study. Two flow conditions are 
examined for each geometry, for a total of four cases. 
In addition, some Mach 6 conditions were revisited. 
The freestream Reynolds numbers for each case, as 
well as dimensional freestream conditions, are listed 
in Table 3. Since no data are available regarding the 
freestream disturbance conditions in the Langley 31” 
Mach 10 facility, a freestream turbulence intensity of 
0.9% was assumed for all four of the Mach 10 cases. 



a) Forebody and shoulder 



b) Base plane and sting 

Fig. 8 Influence of grid refinement for case 2- 
Mach 10, Reo = 2.8 x 10 5 

Additionally, the wall temperature in each solution is 
fixed at 300 K. Table 3 also lists the dimensions of the 
two grids used for each condition. 

Grid Refinement Study 

We shall begin by examining the influence of grid 
resolution on the heating predictions for each of the 
four cases. In each case, the results shown are obtained 
using the newly developed instability mechanism de- 
scribed above. Fig. 7 shows heating results for the 
first of the four cases, with separate graphs for the 
forebody /shoulder heating distribution and the base 
plane/sting heating distribution. The results are sep- 
arated in this manner because of the large differences 
in forebody and base plane heating magnitudes. As 
shown in Fig. 7(a), both the coarse- grid and fine-grid 
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.leat-transfer predictions agree t'niriv well with the ex- 
perimental data. This plot also demonstrates that 
the present forebodv results agree wed with the lami- 
nar solution performed bv Horvath and Hannemann." 
The grid employed in Ref. 6 was 771 x 301. Hence, it 
appears that the current forebodv result obtained us- 
ing grid 2 is weil-resoived. Furthermore, the forebodv 
iieating results support the notion that transition to 
turbulence does not take place on the forebodv. 

Fig. 7(b) compares the present results on coarse and 
fine grids along the back plane and sting. A laminar - 
fiow heating result for this condition 1 * is also shown 
in this plot. This result demonstrates the extent to 
which laminar calculations underpredict the measured 
heat flux along the sting. Such behavior is typical of 
all the laminar solutions presented in References 5, 6, 
and 7. Therefore, comparisons with laminar-flow pre- 
dictions along the base and sting will not be shown in 
the remaining figures. 

It is clear that increased grid refinement has a dis- 
tinct effect on the sting heat-transfer distributions 
obtained using the present approach. This trend is 
similar to that previously observed in laminar compu- 
tations for these flows. 5-7 While the magnitude of the 
peak heat flux does not change greatly, the location of 
the peak heating is visibly shifted downstream. This 
movement arises because the fine-grid solution predicts 
a slightly larger recirculation zone: hence, the reat- 
tachment point and location of shear-layer impinge- 
ment are moved downstream. A further consequence 
of this behavior is that the predicted onset of transition 
(denoted by the vertical arrows along the horizontal 
axis) is also moved further downstream. This plot also 
shows that increased grid refinement tends to improve 
the agreement between the computation and experi- 
mental measurements for this case. 

A comparison of the coarse-gri d an d fine-grid results 
for the second case is shown in Fig. 8. As in the previ- 
ous case, the forebodv heating distribution compares 
well with experimental data and previous laminar com- 
putations. The sting heating results (Fig. 8(b)) show 
similar trends to the first case, in that the refined-grid 
solution predicts a peak heating located farther down 
the sting. In this case, though, the peak heating mag- 
nitude is not as well-predicted as in the first case; the 
discrepancy between the fine-grid heat flux and the 
experimental data is still roughly 30%. It is not im- 
mediately clear why the agreement for this case is not 
as good as that seen in the other solutions. 

For the third case, Fig. 9 shows that the compar- 
isons are fairly good for both the forebody and sting 
heating predictions on the fine grid, with the predicted 
forebody results within or slightly outside of the ex- 
perimental error bounds. Note that the grid used by 
Hollis and Perkins 7 for cases 3 and 4 was a 125 x 357 
grid. The good agreement between the fine-grid re- 
sults obtained here and the laminar solutions of Hollis 



a) Forebody and shoulder 



b) Base plane and sting 

Fig. 9 Influence of grid refinement for case 3- 
Mach 10, Rep = 3.15 x 10 6 

and Perkins indicates grid convergence on the forebody 
for these cases, as well. Along the sting (Fig. 9(b)), 
we again see that increased grid refinement yields a 
later heating peak, as well as delayed onset of transi- 
tion. The refined-grid result also predicts a slightly 
higher peak heat-transfer rate. As a consequence, 
most of the features of the sting heating distribution 
are quite well-predicted, including the initial rise in 
heat flux corresponding to shear-layer impingement. 
The anomalous "kinks” in the heating distributions for 
this and the next case are consequences of a discontin- 
uous slope in the streamwise gridlines, and should not 
be confused with convergence problems. 

The trends in case 4 are similar to those observed 
in the previous cases. As Fig. 10(a) demonstrates, 
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a) Forebody and shoulder 



b) Base plane and sting 

Fig. 10 Influence of grid refinement for case 4— 
Mach 10, Re© = 8.23 x 10 4 


agreement on the forebody is again quite good. Agree- 
ment along the sting is also good; the fine-grid solution 
shown in Fig. 10(b) compares well, in terms of both 
magnitude and location, with the measured peak heat- 
ing rate along the sting support. Once again, the 
coarse-grid solution predicts the peak heating levels 
fairly well, but produces a peak heating location visi- 
bly upstream of that seen in the data. 

Based on the trends seen in these results, as well 
as the trends in other researchers 1 laminar computa- 
tions, further grid refinement would most likely push 
the peak-heating locations slightly further back on the 
sting without greatly altering the magnitudes of those 
peaks. 


Comparison to Original Transition Model 

The next two figures compare representative coarse- 
grid results obtained using the new instability mech- 
anism to solutions using the earlier approach. While 
the results previously presented clearly illustrate the 
impact of grid refinement on the wake solution, com- 
parison of coarse-grid results is still useful, since grid 
refinement does not appear to yield very large changes 
in the magnitudes of the predicted peak heating rates. 

Fig. 1 1 compares sting heating results obtained us- 
ing the old and new instability mechanisms. It is clear 
from this graph that the new mechanism provides a 
greater heating augmentation, and hence better agree- 
ment with the data. 

Comparisons between the old and new approaches 
for the third case are shown in Fig. 12. It should 
be pointed out that the solution using the earlier in- 
stability mechanism was actually obtained on a grid 
with greater normal resolution. In light of the grid- 
refinement study results presented earlier, though, 
comparisons between these two results may still be 
useful. The plot demonstrates again that the new 
mechanism provides substantially greater peak heating 
along the sting support, resulting in improved agree- 
ment with the experimental data. 

While these results indicate that the new instability 
model performs better for the Mach 10 cases exam- 
ined in this work, it is reasonable to ask whether it 
also performs well for other Mach numbers. The new 
mechanism was calibrated using the data from the 
earlier Mach 6 case at a freestream Reynolds num- 
ber of 2 x 10 6 , so good agreement for this condition 
should be expected. However, comparison of the mod- 
els at other conditions is certainly worthwhile. The 
next plot, Fig. 13, compares the two mechanisms for 
another of the Mach 6 flowfields of Horvath et al. 5 
at a freestream Reynolds number of approximately 
0.5 x 10 6 (dimensional conditions listed in Table 2). 
Again, the 189 x 75 grid is used for both solutions. 
This graph demonstrates that the new model yields 
improved heating predictions for this case as well. 

The results shown in Figures 11 through 13 are 
believed to be sufficient to demonstrate the improve- 
ments obtained with the new mechanism. It is inter- 
esting to note that, in all but one of these comparisons, 
the predicted location of transition onset (which cor- 
responds to the location of minimum skin friction) was 
the same for both models. Increases in grid refinement 
clearly have a much stronger impact on the location of 
onset than details of the instability mechanism for the 
cases considered here. 

Concluding Remarks 

The combination of transition and turbulence mod- 
els used in this work was shown to provide improved 
agreement between the numerical heat-transfer predic- 
tions and available experimental data for the Mach 
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Fig. 11 Comparison of earlier and current 
approaches for case 1-Mach 10, Reo = 1 x 10 6 
(coarse- grid solutions) 



Fig. 12 Comparison of earlier and current ap- 
proaches for case 3 -Mach 10, Reo = 3.15 x 10* 
(coarse- grid solutions) 

10 cases considered. In general, the present method 
yielded good comparisons to the experimental sting 
heating data for a wide range of freestream Reynolds 
numbers. 

Coarse-grid solutions produced good predictions of 
the magnitude of the peak sting heating rates, al- 
though the locations of these peaks were consistently 
predicted upstream of the peaks observed in the mea- 
surements. Increasing grid resolution, however, tended 
to improve the agreement even more by shifting the 
predicted peak locations further downstream. 

The overall agreement between computation and 
data for the transitional-flow approach described 
herein is much better than that seen in laminar com- 

10 



Fig. 13 Comparison of earlier and current ap- 
proaches for Mach 6 case, Reo = 0.5 x 10® (coarse- 
grid solutions) 

putations for the same flow conditions. This improve- 
ment lends further support to the notion that transi- 
tion to turbulence does occur in the blunt-body flows 
examined in this study. 

The results presented here indicate that treating 
transition in a turbulence-like manner is a viable ap- 
proach for studying transitional/ turbulent high-speed 
flows. Moreover, it illustrates the importance of em- 
ploying the correct transition mechanism in formulat- 
ing the closure model. Thus, although stability codes 
did not play a direct role in the present calculations, 
linear stability theory results were very important in 
developing the present instability mechanism. 
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An examination of hypersonic turbulent flows where shock/boundarv-layer interaction takes place is 
undertaken. A k-C stress model is developed and used to consider two experiments: a Mach 5 
shock/boundary -layer interaction and a Mach 9.22 compression corner. The stress model, like the two- 
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requires the simultaneous consideration of additional equations governing the variances of density and 
temperature as well as their dissipation rates. 


Nomenclature 

C Skin-friction coefficient 

k Turbulent kinetic energy 

Pfc Production of turbulent kinetic energy 

p Pressure 

Q Reynolds- or time-averaged value of Q 

Q Favre-averaged value of Q 

q Heat-transfer rate 

Re Reynolds number 

T Temperature 

t- Laminar stress tensor 

U Velocity magnitude 

n ; Velocity vector 

3 Flow deflection angle, degrees 

o. ■ Kronecker delta 

:: Turbulent dissipation 

Dynamic viscosity 
v Kinematic viscosity 

p Density 

T t ; Reynolds stress tensor 

in- Vorticitv vector 

^ Enstrophy 

Subscripts 

d Deviatonc component 

T Turbulent value 

u 1 Wall value 

Freestream value 
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Introduction 

The interaction between a shock wave and turbulent 
boundary layer is one of the most significant fluid-dynamic 
problems currently facing designers of high-speed vehicles. 
Shock/boundary-Iaver interactions may be found in 
airbreathing inlets, at wing-bodv junctures, and on 
deflected control surfaces. Accurate characterization of 
these interactions, including resultant heat-transfer rates 
and prediction of separated regions, is very important for 
the design of efficient hypersonic flight configurations. 

Attempts at modeling these types of interactions in 
hypersonic flows have not yet met with the success that has 
been achieved at lower speeds. The object of this 
investigation is to examine some of the limitations of 
existing model approaches, dispel some current 
misconceptions, and suggest ways to develop the models 
necessary for the accurate description of hypersonic 
turbulent flows. This discussion will be carried out in 
coniunction with two experiments that involve shock- 
induced separation. The first 1 is a well-documented 
experiment by Schiilein et al. involving the interaction of 
an oblique shock wave with a Mach 5 turbulent flat-plate 
boundary layer. The second is an older experiment by 
Coleman et al. 2 at Mach 9.22. This particular dataset 
survived the scrutiny of Settles and Dodson. 2 The first 
experiment was the subject of a recent investigation by 
Lindblad et al. 4 using a variety of two-equation and 
algebraic stress models. The second experiment was 
examined by Grasso and FalconF using a k-z model. 

The k-C model 6 will be used as the basis of this study. 
The model performs well for both free and wall-bounded 
shear flows with and without separation, It is free of 
damping and wall functions and is coordinate-independent. 
Moreover, modeled correlations are tensorially consistent 
and invariant under a Galilean transformation. The model 
is limited :o flows where Morkovin's hypothesis 7 is 
expected to hold. 

ft has been suggested 8 * 9 that the inability of existing two- 
eauation models to predict surface properties is a result of 
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the consideration that the turbuient Kinetic energy 
production. P (1 is inversely proporuonui to me end spacing 
normal to the shock. When this behavior is coupled with 
the tact that two-equation models do not accurately predict 
normal stresses even at low speeds, the production term in 
the k equation can assume inaccurate values. Incorrect 
prediction of this term will invariably lead to incorrect 
predictions of surface properties. \s a result of this 
behavior, it is now a common practice to limit i\ to some 
multiple of the dissipation rate. That is. the definition 


is often replaced with 


Modeling Approach 

In this work, we solve the Favre-averaged full Navier- 
Stokes equations for compressible two-dimensional flows. 
Closure is obtained througn use of the two-equation k 
turbulence model and one of two methods for determining 
the normal Reynolds stresses. 

k-t Model Equations 

Turbulence length and time scales are provided in this 
work by the compressible model first proposed by 
Robinson. Harris and Hassan. 11 The final version of the 
model equations is given in Ref. 6 and is also listed below. 
Model constants are listed in Table 1. The exact 
unmodeled k-t equations are given in the Appendix. 
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where ~ :j is the Reynolds stress tensor unu a ;> a user 
specified constant. 

When a stress model is implemented, it is important to 
note that the integrated production of turbulent kinetic 
energy across a normal shock can be written as 
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w'hich is finite and gnd-independent. This result would 
suggest that the key to accurately predicting surface 
properties in shock-dommated flows is to use stress models 
rather than two-equation models. 

Based on the above considerations, a k-t based stress 
model is developed and used to study the experiments of 
References 1 and 2. The stress model was tested to ensure 
that it reproduces data for the incompressible flow over a 
fiat plate and in its wake. Contrary to the prevailing 
wisdom, it is shown that the two-equation model behaves 
in a manner similar to that of the stress model: both 

perform well for the Mach 5 experiments, and both did 
poorly tor the Mach 9.22 experiments. This finding was 
not much different from that of Zha and Knight. 10 who 
used a stress model to study a crossing-shock problem: its 
prediction of surface properties was not much better than 
the underlying two-equation model. 

Based on this investigation, it appears that accurate 
predictions of strong shock/boundary-layer interactions 
requires the relaxation of Morkovins hypothesis. This 
relaxation, in turn, requires supplementing existing 
turbulence models with equations governing the variances 
ot density and/or temperature together with their respective 
dissipation rates. 
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As indicated above. Morkovins inpothesis is employed 
n deriving the model equations 1 4) ana ifi from me exact 
equations governing k and 4 ;Eqs. A. I and A. 2 in the 
Appendix). According to this hypothesis, the pressure and 
total-temperature fluctuations are negligible tor non- 
hvpersonic boundary' layers with conventional rates of heat 
transfer. Thus, all correlations involving in Eqs. (A.l) 
and (A. 2) were ignored when deriving above model 
equations. Traditionally. Mach (A/1 5 has been chosen as 
the Mach number below which Morkovins approximation 
is expected to hold when considering wall-bounded shear 
Hows. Thus, based on the above consideration, one should 
expect that current k-C formulation in either a stress or two- 
equation version, to perform weil for wail bounded flows at 
A/=5 or less, but not for higher Mach numbers. Above 
A/=5, amplitudes of density and temperature fluctuations 
become significant and can not be ignored. Therefore, at 
hypersonic Mach numbers, equations governing variances 
of density and temperature together with their dissipation 
rates must be part of a turbulence model. 

The stress model developed here again invokes 
Morkovin’s approximation. This will facilitate assessment 
of penormance of stress and two-equation models when 
used to solve a problem where Morkovin s approximation 
is expected to hold and where it is not expected to hold. 

Reynolds Stress Closure 

The Boussinesq approximation, i.e.. 

— f 1 2 

T.. 3 - p u’u’j = u r , 2s ij - j 5 ,y.v u - -=• 0 ,, pA i 6) 

lias been widely used in turbulence modeling. This 

assumption generally yields good results for the snearing 
stress but not for the normal stresses, as is indicated in Fig. 

1. This plot compares predicted root-mean-square (RMS') 
velocity fluctuations using the k-C model to the 
experimental data obtained bv Klebanorr ,: for 

incompressible flow over a flat plate. Application of the 
thin-shear-layer approximation to Boussinesq’s assumption 
yields purely isotropic normal stresses, such that 

Tu: ~ ^ vv ~ ^ :z ~ 7 P^ 
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Consequently, all three RMS fluctuation intensities follow 
a single curve, as demonstrated in Fig. 1. Cleam . the two- 
equation approach yields inaccurate results for the normai 
stresses even for this very simple flowtield. An alternative 
approach is to pursue a ; tress formulation. A differencial 
stress model adds five equations to a two-equation model. 
This additional complexity and cost has discouraged 
potential users in the past and shifted attention to Algebraic 
Reynolds Stress Models (ARSM). These algebraic stress 


models arc derived from differential stress models by 
making two important assumptions: diffusion is negligible 
and turbulence i> m equilibrium. Both of these 
assumptions are invalid for separated flows. Therefore, 
one should not expect ARSM to improve our predictive 
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Fig. 1 Comparison of computed and experimental 
velocity fluctuations for flat-plate experiment of 
Klebanoff-two-equation model 

capabilities of surface properties for such flows. This 
observation was confirmed by, among others, a recent 
investigation conducted by Rizzetta 13 in which he studied 
shock wave/boundary layer interactions using three ARSM, 
two two-equation models and an algebraic model. He 
concluded that ARSM offered little improvement over the 
other models. 

A stress model based on the k-C formulation is developed 
here. In order to avoid the expense of employing a full 
stress model, all calculations presented here were based on 
a Partial Differential Stress ( PDS) model. In this approach, 
the shear stress is calculated from a Boussinesq 
approximation while the normal stresses are calculated 
using modeled transport equations for these quantities. 
Thus, whether one deals with two- or three-dimensional 
flows, the overhead for this model is two additional field 
equations. 

To formulate this new approach, we begin with the exact 
Reynolds stress equation 14 
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The second-order tensors e,y and D.. may be decomposed 
into isotropic and deviatoric components. Thus. 
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The exact Reynolds stress equation can then be rewritten as 
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Upon contraction. Eq. (8) reduces to Eq. (A.l) in the 
Appendix. Therefore, assumptions that were used in 
modeling the k equation can be adopted here. This ensures 
that the underlying k formulation is the same for both the 
two-equation and the stress model. An added advantage is 
that the only term that requires modeling is o,y . Hence, the 

model Reynolds stress equation is 


Consistent with the practice employed in developing the 
k-Z model. 0 11 the quantity o ( j will be modeled to be a 

linear combination of the anisotropy tensor. by , and its 
derivatives. The resulting expression can be written as 
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Note that each model term contracts to zero independently. 

The constants C L - C 4 can depend on the invariants of 
b { j and P k /pe. As a first step they will be treated as 

constants. They will be determined from rapid distortion 
theory, which gives C 2 , and the rest from Laufer’s 

experiment in channels. 15 Rapid distortion theory gives 
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Laufer’s measurements give, in the log-law region where 
Pk/ P £ 38 j * 

b { j = 0.22. bi2 = “ 0.15. d 12 =-0.15 (13) 

Equations 1 12) and (13) yield 

C I =2.513. C 3 =1.473. C 4 =1.07 (14) 

An adjustment of the above model constants is required 
for situations where P k / pe * l . Abid and Speziale 16 have 
suggested that second order closures should yield 
(10) approximately equal equilibrium values for by for both 
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Homogeneous shear rlows and ihe log-jaw region. This 
goai can be achievea by a simple moditlcation or' c\ to 

C { = 2.0 + 0.5 13 (P; / p£ ) ( 15 ) 

Note that Cj here is 2Cj in Ret. 16.) The above choice 
replaces Rotta's constant hy 2.0 and results in 
homogeneous shear flow values of b :j that are given by 

Eq. \ 13) for all values of P k / p£ ( see Eq. ( 1 1) in Ref. 16.] 
The constant C 3 was chosen so as to reproduce correct 
skin friction for a flat plate. Thus, the model constants for 
low speed flows in the absence of adverse pressure 
gradients can be summarized as 

C[ = 2.0 + 0.5I3(/ype>. C : =0.8. C 3 = l .473. 
C 4 =1.07, C 5 =0.1 (16) 


Again, consistent with the development of the k-C model, 
the role of compressibility and adverse pressure gradients 
had to be addressed. From study of shock wave/boundarv 
layer interaction. C\ was multiplied by the factor 



where 



The model constants (3 p and P p . as well as the other 
stress model constants, are listed in Table 2. 


Table 1 k-C / PDS model constants 


Constant 

Value 

C M 

0.09 

a% 

0.35 

Pa 

0.42/0.39 

05 

2.37 

06 

0.1 0/0. IS 

07 

1.50 


0.065 

O r 

0.07 

o p 

65.0 

I/O;. 

1.80 

I/O; 

1.46 

<5 

0.10 

Ci 

0.60 

c k 

1.00 
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r i 


Table 2 PDS model constants 


Constant 

Value 

Cl 

2 + 0.513 P k /(pe) 

C i 

0.8 

C 3 

1.473 

C 4 

1.07 

C 5 

0.10 

0, 

0.18 

0P 

0.57 


Note that, in developing this model, we preserved 
all the desirable features of the two-equation model, 
including avoidance of wall or damping functions and 
coordinate independence. Moreover, all correlations are 
tensoriallv consistent and Galilean invariant Further, as 
shown in Table 1 the majority of the k constants remain 
unchanged for the PDS model. The constants f } 4 and P$, 
however, are modified slightly as shown. 

Numerical Method 

The modeling proposed here has been incorporated into 
Olynick and Hassan's 17 two-dimensional/axisymmetric 
implicit solver for hypersonic flows. This algorithm solves 
the governing equations for 5 -species air in 
thermochemical nonequilibrium. The solver uses Roe's 
flux -difference splitting 18 for the in viscid flux, extended to 
higher order using MUSCL variable extrapolation 19 with a 
minmod slope limiter. Time integration is accomplished 
using the diagonal implicit variant of Yoon and Jameson's 
Lower-Upper Symmetric Gauss-Seidel method, 20 which 
only requires the inversion of diagonal matrices. This 
property is attractive for nonequilibrium flows, where a 
large number of partial defferential equations must be 
solved. 

Since the current flowfields under investigation are 
perfect-gas. vibrational relaxation and chemical reactions 
are disabled to obtain the results in this study. 
Additionally, the high-temperature transport-property 
calculations originally used in the code were replaced with 
Sutherland's law for viscosity and a constant laminar 
Prandtl number of 0.72. Closure for the Favre-averaged 
energy equation is accomplished using a constant turbulent 
Prandtl number of 0.89. 

Results and Discussion 

Before considering high Mach number flows, it is 
worthwhile to examine the performance at the PDS model 
for low speed flows in the absence of adverse pressure 
gradients. This study will help evaluate the choice of 
model constants indicated in Eq. (16). A flat plate and its 
wake were selected for this validation. Of the many results 
that were generated, two critical comparisons will be 
presented. Figure 2 compares computed RMS velocity 
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iluctuations with Klebanorfs experiment.* 2 while Fig. 3 
compares computed wake growth rate with those of the 
two-equation model 0 and the experiments of Pot :i and 
Wevgandt and Mehta. :: In Fig. 3, b is the half width. 0 is 
the momentum thickness and .v is the distance along the 
wake. As is seen in the figures, good agreement is 
indicated in both cases. 



Fig. 2 Comparison of computed and experimental 
velocity fluctuations for flat-plate experiment of 


Klebanoff-PDS model 


700 h 



600 f- 


available data sets for code validation. However, the 
experiment met all their guidelines. The work of Coleman 
et ai. was on their list of recommended experiments. Note 
that, for each of the computations presented here, the 
leading edge of the flat plate or compression ramp is not 
included in the Navier-Stokes computations. Instead, 
inflow profiles are obtained separately through boundary- 
layer solutions for the corresponding zsro-pressure- 
gradient flat-plate flow. 

Computations for Schiilein et ai. Cases 

First, we shall examine a Mach 5 shock/boundary-layer 
interaction modeled using the experimental setup shown in 
Fig. 4. A shock generator is mounted on the upper tunnel 
wail, and the resultant oblique shock wave interacts with 
the turbulent boundary layer growing on the flat plate along 
the lower tunnel wall. Depending on the magnitude of the 


SHOCK GENERATOR 



Fig. 4 Schematic of shock/boundary-layer interaction 
experiment performed by Schiilein et ai. 
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Fig. 3 Comparison of k-C, and PDS results for 
incompressible flat-plate wake 

We shall consider next two different classes of 
shock/boundary-layer irteraction flowfields. including 
Mach 5 impingement of an oblique shock wave and the 
flow past a compression comer at Mach 9.22. These 
computations correspond to experiments conducted by 

Schiilein et al.^ and Coleman et at..- respectively. It 
should be noted that the first experiment was completed 
after Settles and Dodson completed their examination of 


deflection angle (3, flow separation may or may not occur at 
the location of shock impingement. We consider two 
values of 3 which result in separated flow. The oblique 
shock is enforced numerically by enforcing velocity, 
pressure, temperature and density boundary conditions 
along the top of the computational domain that correspond 
to the conditions behind the oblique shock. It was 
indicated in Ref. 4 that an 80 x 80 grid resulted in a well 
resolved solution. Therefore, in order to ensure a grid- 
resolved result, all present computations were conducted on 
a 141 x 141 Cartesian grid w ith constant spacing in the x 
direction and geometric spacing in the v direction. 

Before detailed comparisons are undertaken, the 
turbulence model of Ref. 4 will be reviewed. Lindblad et 
al. 4 employed a number of two-equation, ARSM and 
differential stress models in their investigation. We will be 
comparing here with what they considered to be their best 
performing model, which happened to be an Explicit 
Algebraic Reynolds Stress Model (EARSM) with the 
underlying two-equation model being the k- o model of 
Wilcox. This EARSM is based on the ARSM of Rodi 23 
and the pressure-strain model of Launder et al. 24 A near 
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wail treatment, which included damping functons. was then 
added- 5 *- 6 to allow the Reynolds stress anisotropy to 
behave in accordance with experiment and direct numerical 
simulation. Evidently, 4 compressibility corrections as 
suggested in Refs. 2~ and 28 were deemed unimportant and 
thus, none were included. 

Fig. 5 compares nondimenstonai surface pressure and 
skin-friction results ootained with the standard k-C model 
and the PDS model to the experimental measurements of 
Schulein et al. for (3 = 10°. Also shown are the results 
obtained by Lindblad et al. 4 using an EARSM for the same 
case. Fig. 5(a) shows that the present computations both 
predict the downstream pressure quite well, but offer 
differing results in the interaction region. These 
differences appear to arise because of a considerably 
smaller separation zone predicted by the PDS model, as 
indicated in Fig. 5(b). Outside the interaction region. 



a)Pressure distribution 



b) Skin-friction distribution 


Fig. 5 Results for experiments of Schulein et al.< (3 = 10° 


however, all models compare quite well to the 
experimental skin-friction data. For this case, the 
agreement appears to be better than that provided by the 
EARSM computation. 

Fig. 6 shows that the agreement is generally poorer for 
the case where (3 = 14°. For instance, the present two- 
equation computation predicts the general behavior of the 
pressure distribution quite well: however, it underpredicts 
the pressure in the interaction region and shows an 
overprediction in the downstream equilibrium region. The 
PDS model better predicts the downstream pressure but 
does not agree well in the interaction region. For the skin 
friction, the PDS model appears to match the skin-friction 
data better than the standard k-C approach, but predicts later 
separation and reattachment. It is also clear that both the 
standard k-C model and the EARSM calculation yield large 
overpredictions of the skin friction in the recovery 



a) Pressure distribution 



b) Skin-friction distribution 


Fig. 6 Results for experiments of Schulein et aL, 3 = 14° 



region. Overall, the k-Z and PDS calculations provide 
similar degrees of predictive accuracy for these Mach 5 
cases. 

It is to be noted that, in ail of our calculations presented 
here, no attempt was made to limit the ratio of production 
to dissipation. Detailed investigation of resuits from the 
two-equation model indicated that, in regions where 
production over dissipation was unsuallv high, production 
was dominated by normal stresses. This may be the reason 
why skin friction was overpredicted in the recovery region 
by the two-equation model. 

Computations for Coleman et al. Cases 

The second set of experiments considers the tlow past a 
two-dimensional compression comer in a Mach 9.22 
freestream. For these cases, we shall compare 
computational results to experimental surface-pressure and 
heat flux results obtained by Coleman et al. 2 A 
representative computational grid for one of the two 
deflection angles considered is shown in Fig. 7. The 

minimum grid spacing at the wail decreases linearly from 
the inflow plane to the comer, after which point it is 
constant. Two flow deflection angles. 3 = 34° and 3 = 38°, 
are examined for this case as well; flow separation is 
expected for both deflection angles. Both the 34- and 38- 
degree grids employed 241 points in the streamwise 
direction and 141 points in the normal direction. This is 
comparable to 256x128 grid employed in Ref. 5. which 
resulted in a grid-independent solution. 



Fig. 7 241x141 grid for Coleman et al. experiment. 
3 = 34°. Every other grid point omitted for clarity. 

The calculations of Ref. 5 were based on a two-equation 
k-E model with a number of compressibility corrections 
involving terms appearing in the k equation, and the 
Karman constant. Moreover, the eddy viscosity definition 
was adjusted so as to limit the turbulent length scale. 

As indicated earlier, because all correlations involving 
p were set to zero, our model is not expected to perform 
well for the Af=9.22 case. The need to include temperature 


and density fluctuations can be seen when one considers, as 
an example, the pressure dilatation term pu\ j . Since 

p = p +■ p = Rp{ T ^ T') 

p = /?(pT - pn 

and 

p'u'u = R[T p'u'i i + p T'u ' u ] , (19) 

it is not obvious why the above term should be modeled in 
terms of turbulent kinetic energy production P k and its 

dissipation rate pe , as suggested by Sarkar 28 and used in 

Ref. 5. We believe that appropriate modeling of terms 
appearing in Eq. (19) must depend on the variances of 
temperature and density and their dissipation rates. 

We shall first compare predictions for inflow properties 
with these of Ref. 5 and experiment. These results are 
shown in Fig. 8. Fig. 8(a) shows that both the ifc-£ solution 
and the PDS solution match experimental velocity-profile 
data better than the earlier calculation. In the case of the 
Mach-number distribution (Fig. 8(b)), the Grasso and 
Falconi result agrees better with the data in the outer 
region, but the present computations are superior in the 
near-wall region. 

It is seen from Fig. 8 that, in general, current predictions 
are better than those of Ref. 5 and in good agreement with 
experiment. The question arises as to why a theory that is 
not supposed to be valid in this Mach number range gives 
good agreement with experiment. This is because 
Morkovin’s hypothesis requires conventional heat transfer 
rates which is the case for the inflow region. This is to be 
contrasted to the comer region, where high Mach numbers 
and unconventional rates of heat transfer exist. 

The next two figures show pressure and heating results 
for the two deflection angles considered here. Fig. 9 
compares the present results to the computations of Grasso 
and Falconi and experimental data for 3 = 34°. As 
demonstrated in Fig. 9(a), the present k-Z> computation 
appears to provide the best comparison with the 
experimental pressure data for this case, whereas the PDS 
result provides the worst agreement. On the other hand, the 
k-E solution by Grasso and Falconi agrees fairly well with 
the experimental heat-transfer data, but both present 
computations provide substantial overpredictions of the 
peak heat-transfer rate. Note that the k-C result agrees with 
the data fairly well upstream and downstream of the 
interaction region; however, the predicted peak heat flux is 
more than twice that observed in the experiment. The PDS 
heating result overpredicts the heating in and downstream 
of the interaction region. 

Similar behavior is observed for 3 = 38° in Fig. 10(a); 
both the k-Z and k-E solutions predict the peak pressure 
magnitude quite well. Once again, however, the 
computation by Grasso and Falconi provide reasonable 
predictions of the heat-transfer rate, while both present 
computations overpredict the peak heating substantially. 
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Moreover, the PDS result again agrees poorly with the data 
downstream of the interaction region. 

Concluding Remarks 

This work extended the two-equation k-C turbulence 
closure model to a stress model and, in the process, 
preserved all the desirable features of the original model, 
i.e., no wall or damping functions, coordinate 
independence, Galilean invariance and censorial 
consistency. 

The above results illustrate the importance of using a 
turbulence model that reflects the correct physics of the 
problem. The k -£ model worked for the Mach 5 shock 
wave/boundary-layer interactions and for the 9.22 flat-plate 
boundary-layer computations. However, it did not provide 
accurate predictions for the separation zone in the Mach 



a) Velocity profile 



c) Mach-number profile 


Fig. 8 Comparison of inflow profiles for experiment of 
Coleman et al. 


9.22 compression-comer flowfields. As a result, we 
believe that the key to developing accurate models for 
hypersonic turbulence is to relax Morkovin's 
approximation and provide adequate modeling for terms 
neglected in Eqs. (A.l) and (A.2). Replacing a two- 
equation model with a stress model is not going to lead to 
improved predictions if the underlying physics is 
inaccurate. 

Another conclusion that may be drawn from this work is 
that, if the objective is to calculate surface properties, then 
two-equation models perform as well as stress models as 
long as the production term P k is not dominated by 

contributions of normal stresses. 

Finally, the fact that the PDS model underperformed the 
two-equation mode! for the Mach 9.22 compression ramp 
was unexpected. Future work should include comparison 



a) Pressure distribution 



b) Heat-transfer distribution 


Fig. 9 Results for experiments of Coleman et aL, 0 = 34° 
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oi the PDS model with a full stress mouei. 
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Appendix: Exact Equations Governing k and £ 

The exact Favre averaged compressible turbulence 
kinetic energy and enstrophy equations can be written as 
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